/*	Interactive Gravitational Simulator
	Project: IGS Core
	File: Plummer.cpp
	Author: Mike Bantegui <mbante2@gmail.com>, Hofstra University 
	Copyright (C) 2012 - 2013 Mike Bantegui

	This file is part of the IGS software package for simulating
	the N-Body problem in real-time. This was developed as part
	of a departmental honors thesis project at Hofstra University
	in the Spring 2012 semester.

	IGS is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation; either version 2 of the License, or
	(at your option) any later version.

	IGS is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/

#define _USE_MATH_DEFINES
#include <cmath>
#include "Body.h"
#include "Plummer.h"

namespace IGS
{

// Destructs the PlummerModel
PlummerModel::~PlummerModel()
{
}

// Creates a system of N bodies in a cold collapse unit sphere configuration
void PlummerModel::Fill(Body *bodies, const int N)
{
	const double mass = 1.0 / N;
	const double scale = 16 / (3 * M_PI);
	double cumulativeMin = 0;
	double cumulativeMax = mass;

	// See implementation at
	// http://www.artcompsci.org/kali/vol/plummer/ch02.html#rdocsect9
	for (int i = 0; i < N; i++)
	{
		const double cumulativeMass = Random(cumulativeMin, cumulativeMax);
		cumulativeMin = cumulativeMax;
		cumulativeMax += mass;

		double x, y;

		do
		{
			x = Random(0, 1);
			y = Random(0, 0.1);
		} while (y > x * x * pow(1 - x * x, 3.5));

		const double r = 1 / sqrt(pow(cumulativeMass, -2.0 / 3.0) - 1.0);
		const double v = x * sqrt(2 * pow(1 + r * r, -0.25));

		const Vector position = RandomSpherical(r) / scale;
		const Vector velocity = RandomSpherical(v) * sqrt(scale);
		bodies[i] = Body(mass, position, velocity);
	}
}

}